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ABSTRACT 

We describe the time-dependent radiation transfer in blazar jets, within the 
internal shock model. We assume that the central engine, which consists of a 
black hole and an accretion disk, spews out relativistic shells of plasma with dif- 
ferent velocity, mass, and energy. We consider a single inelastic collision between 
a faster (inner) and a slower (outer) moving shell. We study the dynamics of the 
collision and evaluate the subsequent emission of radiation via the synchrotron 
and synchrotron self Compton (SSC) processes after the interaction between the 
two shells has begun. The collision results in the formation of a forward shock 
(FS) and a reverse shock (RS) that convert the ordered bulk kinetic energy of the 
shells into magnetic field energy and accelerate the particles, which then radiate. 
We assume a cylindrical geometry for the emission region of the jet. We treat the 
self-consistent radiative transfer by taking into account the inhomogeneity in the 
photon density throughout the region. In this paper, we focus on understanding 
the effects of varying relevant input parameters on the simulated spectral energy 
distribution (SED) and spectral variability patterns. 

Subject headings: BL Lacertae objects: general — Galaxies: jets — Hydrodynam- 
ics — Radiation mechanism: non-thermal — Radiative transfer — Relativistic 



processes 
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Introduction 



Blazars are a class of extreme active galactic nuclei (AGNs). They consist of flat- 
spectrum radio quasars (FSRQs) and BL Lac objects, and are known for their rapid vari- 
ability, at all wavelengths, on timescales of months to a few days and e ven to less than 
an hour in some cases (e.g., Gaidos et al. 1996: Catanese & Weekes 1999; Aharonian et al. 



20071 ; Albert et al.l 120071 ) . Such variability is indicative of highly violent phenomena and the 
presence of ultra-relativistic particles in the jets of blazars. The SED of blazars consists of 
two broad spectral bumps that are associated with synchrotron radiation for the low-energy 
component and inverse Comptonization for the high-energy one. The spectral variability 
patterns and SEDs provide valuable information regarding the acceleration of particles and 
the time-dependent interplay of various radiation mechanisms responsible for the observed 
emission. It is, therefore, imperative to develop a theoretical approach that is capable of 
simulating these observables to gain a deeper understanding of the physics of blazar jets. 



Many authors have studied bla z ar jets in the past (see. e.g. jMaraschi. Ghisellini fc Celotti 



1992 



Dermer fc Schlickeiserl Il993l ; IBottcher. Mause fc Schlickeiser 



1997 



Joshi fc Bottchcr 



20071 ) in the framework of one-zone models. In those works, a spherical emission region 
has been considered and a power-law energy distribution of injected electrons and positrons 
assumed. The evolution of the particle population and the subsequent photon emission is 
solved numerically in a time- dependent manner in order to generate theoretical light curves, 
spectral variability patterns, SEDs, and other features to obtain a satisfactory fit to the 
observed data of a given source. Such theoretical efforts have successfully reproduced the 
observed SED and the spectral variability patterns of various blazars. They have provided 
information about the relevant physical parameters of the emission region, but the location 
and the mode of particle acceleration in the jets of such systems still need to be explored in 
more detail. Most of the theoretical approaches have assumed a spherical geometry for the 
emission region, continuous injection of particles, and homogeneous population of both elec- 
trons and photons throughout t he emission region (e.g. lSikora et al.lll994j ; lBloom fc Marscher 
19961 ; IBottcher fc Bloom! 120001 ) . Such assumptions are good only to a first approximation 
and do not completely account for the observations. Thus, a more sophisticated approach is 
required to treat the characteristics of the energized plasma, such as electron energy distri- 
bution and magnetic field, not as free parameters but as products of global parameters. This 
would enable us to link the observed emission properties to the physics of energy transport 
inside a jet. 

In the case of leptonic jet models, various authors have t reated the que s tion o f particle 
acceleration and their subsequent injection in different ways. ISokolov et al.l ( 120041 ) consider 
the formation of a system of standing oblique shock waves in the jet that accelerate parti- 
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cles. The local pressure imbalance between the ambient medium and the jet gives rise to 
standing oblique shock waves, which can terminate in a Mach disk, near the axis, that lies 
perpendicular to the jet flow. As the relativistic flow of plasma passes through the disk, the 
plasma electrons accelerate to highly relativistic energies and radiate away their energy via 
synchrotron and inverse Compton processes to produce the observed emission. 



Spada et al.l (120011 ) (hereafter SOI) and lMimica et al.l (120041 ). consider the internal shock 



model to accelerate the particles inside a jet. In this case, shells of plasma with different 
mass, energy, and velocity are ejected intermittently by the central engine. The faster moving 
shells catch up with the slower moving ones to result in collisions that produce shocks internal 
to the jet. These shocks accelerate the plasma electrons to ultrarelativistic energies, which 
lose their energy via various radiative processes to produce the observed radiation. 

All of the above-mentioned approaches consider the formation of reverse and forward 
shocks that lead to particle acceleration in a region inside the jet. The injection of particles 
into the emission region continues until the shocks cross the entire region. The population 
of electrons and photons has been generally approximated to be homogeneous in density 
throughout the emission region, while the geometry of the reg i on ha s been assumed to be 
cylindrical (cubic geometry in case of IChiaberge &: Ghisellinil (ll999[ )). as appropriate for 
internal shocks. 



Graff et al.l (120081 ) invoke a standing or propagating shock in a collimated jet to accel- 
erate particles. They adopt a multi-zone pipe geometry for the emission region to simulate 
the resultant radiation and variability of blazars via synchrotron and SSC processes. Their 
model considers the inhomogeneity in the particle and photon population only in the axial 
direction when calculating the non-local, time-delayed SSC emission of the source. Their 
model does not take into account the volume and angle-averaged photon escape timescale 
for a cylindric al region. Other accele r ation scenarios, s uch as isolated sho cks propagating 



along the jet ( IMarscher &: Gear! Il985t iKirk et al 



1998 



Sikora et al.ll200ll ). oblique shocks 



propagating thr ough a pre-existing je t ( lOstrowski & Bednarzll2002l ). or particle acceleration 
in shear flows (IRieger &: Duffvl 12004 ) have also been considered as plausible acceleration 
mechanisms in blazar jets. 



In a recent approach, iBottcher fc Dermerl (120101 ) employ the internal shock model to in- 
vestigate its effects on the time-dependent spectral evolution and cross-frequency correlation 
of the emission from the jets of blazars. The model considers a one-zone cylindrical emission 
region bounded by forward and reverse shocks that energize the background plasma to very 
high energies. The resultant non-thermal distribution of particles then decays in energy 
through synchrotron, SSC, and external Compton (EC, Compton upscattering of photons 
external to the jet) emission processes. The time-dependent synchrotron and EC spectra are 
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calculated analytically, while the SSC emission is computed semi-analytically, with the cal- 
culations restricted to the Thomson regime. As a result, their model is unable to reproduce 
TeV emission from the high-frequency peaked BL Lac objects (HBLs). 

In this paper, we consider the internal shock scenario to develop a 1-D multi-slice time- 
dependent leptonic jet model for simulating the radiative transfer in relativistic blazar jets. 
We use this scenario to explore the relevant physical parameters of a cylindrical emission 
region that is responsible for the observed radiation. We take into account the inhomogeneity 
in the photon and electron populations throughout the region by dividing it into multiple 
slices and considering time-dependent radiation transfer within and in between each slice. 
The multi-slice scheme with radiation feedback automatically lets us address the non-local, 
time-delayed SSC emission of the source in a self-consistent manner along the length of the 
emission region. It also takes into account the volume and angle-averaged photon escape 
timescale. We calculate the time- dependent synchrotron and SSC emission spectra without 
restricting our calculations to the Thomson regime, hence making it applicable to all classes of 
blazars. We also incorporate the light-travel time delays to correctly calculate the observed 
spectra and light curves in the observer's frame. We do not consider the retarded-time 
photon field individually in the radial direction for calculating the SSC emission neither 
do we evaluate the microphysics of particle acceleration as that is beyond the scope of the 
current work. 

We assume a background thermal plasma in the jet that contains non-relativistic elec- 
trons, positrons, and protons. In §2], we describe the interaction between two shells of plasma, 
responsible for producing internal shocks, and the subsequent treatment of the collision in 
the framework of relativistic hydrodynamics. As the shocks propagate through the plasma, 
the internal kinetic energy of the shocked fluid that is stored in the baryons is transferred 
to electrons and positrons. This transfer of energy accelerates electrons and positrons to 
relativistic and non-thermal energies. In this section, we derive various emission region pa- 
rameters resulting from shock propagation and particle acceleration. In fj3j, we present a 
semi-analytical calculation of the volume and angle-averaged photon escape timescale for a 
cylindrical geometry. In §U we describe the numerical approach that we use for calculat- 
ing various radiative energy loss rates as well as photon emissivities. Since positrons lose 
the same amount of energy as the electrons via the same radiative loss mechanisms, we do 
not distinguish between them throughout the paper. We discuss our multi-slice radiation 
transfer scheme that addresses the inhomogeneity in the density of the photon and electron 
population throughout the emission region, ideals with the relevant light-travel time delays 
that are required to accurately register the radiation in the observer's frame of reference. 
The effects of varying the relevant physical input parameters on the simulated SED and 
lightcurves are discussed in £j6l We present our results of the parameter study in ^3 We 
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summarize our results in §HJ 

Appendices [A] and [B] delineate the details of some of the cumbersome equations used 
in the numerical computation of radiative output. In appendix [UJ we discuss the numerical 
approximation used to calculate the synchrotron photon density in the simulations. 



2. Hydrodynamics of the collision 

In the internal shock scenario, the central engine ejects relativistic shells of plasma with 
disparate mass and energy that travel with different velocities. As a result, faster moving 
shells catch up with slower moving ones, released at earlier times, resulting in inelastic 
collisions. Each collision results in the formation of a FS that propagates into the slower 
moving outer shell and a RS that propagates into the faster moving inner shell (SOI). The 
shocks are separated by a contact discontinuity (CD) across which the pressure of the shocked 
fluids is in equilibrium. As the shocks propagate into their respective shells, they convert the 
ordered bulk kinetic energy of the shocked plasma into magnetic field energy and internal 
energy of electrons. As a result, electrons are accelerated to highly relativistic energies and 
radiate to produce the observed emission throughout the electromagnetic spectrum. 

In the present paper, we consider a single inelastic collision between a slower moving 
outer shell and a faster moving inner shell. We study the dynamics of the collision and 
the subsequent photon emission from non-thermal relativistic electrons after the interaction 
between the two shells has started. The entire treatment of shell collision and shock prop- 
agation is hydrodynamic and relativistic in nature. We do not consider the expansion of 
individual shells or the merged shell along or transverse to the direction of motion. The 
respective exp ansions are neglecte d here because the jet remains well-collimated at sub-pc 



and pc scales (IJorstad et al.l 120051 ) and the collision lasts for a sh ort time during which the 



shell will have a nearly constant width ( IBottcher fc Dermerll2010l ) 



2.1. Collision parameters 

We consider an ejection event of two shells lasting for time t w . The outer shell is ejected 
with rest mass M Q , width A G and a bulk Lorentz factor (BLF) r o . The inner shell is ejected 
with rest mass Mi, width A; and BLF IY The ejection of the two shells is separated by a 
time interval t e . The kinetic luminosity of the entire ejection event is given by L w such that 
the sum of the kinetic energies of the two shells is not more than the total specified through 
L w and t w and satisfies the condition: 
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M o r o c 2 + M^c 2 



L w t w . 



The radii, R, of both shells are taken to be equal to each other and to that of the 
collimated jet. The z-axis is along the axis of the jet and the radius of the jet is perpendicular 
to this axis. At time t — 0, the outer shell is at a distance z Q and the inner shell is at a 
distance z\ from the central engine. In this section, all non-primed quantities refer to the 
rest frame of the AGN (lab frame), all primed quantities refer to the comoving frame of the 
shocked fluid behind the shock front, and all quantities with an ov erline refer to the frame of 
the unshocked fluid. The time of collision ( iKobayashi et al.lll997l ). St can then be calculated 
by 



A ) 



St = ^ - (2) 

and the position at which the two shells would collide, which also defines the location of the 
CD, is given by 



Z\ + cfrSt 



(3) 



Considering the collision of these two shells to be an inelastic one, the BLF, r m , and 
the total internal energy of the merged shell, E- mt , can be obtained by using the conservation 
of momentum and energy (SOI). This yields 



iMjTj + M o r o 

Mi _|_ Ms. 

r t ^ r Q 



(4) 



and 



M iC 2 (r ; - r m ) + m oC 2 (r - r r 



(5) 



The internal energy of both shells before the collision has been assumed to be relatively 
negligible, so that essentially all of the energy is stored as bulk kinetic energy. The efficiency 
of conversion of this energy into the internal energy of the merged shell from a single collision 
is given by 



(Mi + m ) r m 

(MiT; + M o r o ) 



(6) 
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In the comoving frame of the shocked (emission) region, the two shocks appear to move 
in opposite directions, away from the CD, with Lorentz factors r^,^, respectively. In the lab 
frame, the entire shocked material continues to move forward with BLF T^. As a result of 
shock propagation, the emission region keeps increasing until the shocks reach the respective 
boundaries of the merged shell, as shown in Figure [TJ 

The BLF of the entire emission region in the lab frame, r s h, is calculated by using the 
relativistic hydrodynamic Rankine-Hugoniot jump conditions for energy and particle density 
across both shock fronts, and th e condition of equal pressure of the shocked fluids across the 



across Dotn snock ironts, ana tn e condition 01 equal pressure 01 tne snockea nuias across tne 
CD (SOI). The jump equations ( Blandford fc McKeel 1976 ; Panaitescu &: MeszaroslE999l ) for 



the kinetic energy density U' is and the matter density of baryons p' {s in the comoving frame 
of the shocked-fluid frame are given by 



Ul 



fs(rs) 2 



Pfs(rs) 

and 



Pfs(rs) 



* ( r fs(rs) " 1) (7) 



(V fs(rs) + 1 



Po(i) (7-1) 



where p Q ^ = r ~^^a n * s ^ ne ma ^ er baryon density of the preshocked fluid in the comoving 
frame of the unshocked fluid. The quantity 7 is the adiabatic index, whose value is assumed 
to be 4/3 (for an ultra-relativistic plasma). 

The pressure of the shocked fluids is obtained from equations [7] & El and is given by 

Pfe(rs) = ( r fe(rs) - l) (j^rs) + X ) Po(i)C 2 • (9) 

The BLF of each shell relative to the shocked fluids in their respective comoving frames 
can be written as 



r f S (rs) = r (i)r S h (1 - /? (i)Ash) • (10) 

Substituting Eqn. [TD] in Eqn. M and applying the condition of equal pressures yields a 
quartic equation in r sh 



aT* h + &r s 3 h + c r^ h + dT sh + e = 



(11) 
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Fig. 1. — Schematic diagram of the two-shell inelastic collision resulting in a merged shell of 
bulk Lorentz factor (BLF) T m . The forward shock (FS) and reverse shock (RS) fronts start 
to propagate in the outer and inner shells, respectively, creating the forward and reverse 
emission regions. In the comoving frame, the shock fronts move with BLFs T' {s and r^. The 
contact discontinuity (CD) is located at a distance z c from the central engine. The comoving 
pressures p' {s and p' TS of the shocked fluids across the CD are equal. Ultra-relativistic electrons 
are only present in the forward and reverse emission regions. These emission regions keep 
expanding until the shocks reach the respective boundaries of their regions. 
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The values of each of these coefficients are given in Appendix [A] 

As the shocks cross through their respective regions, they compress the plasma present 
in those regions. As a result, the width of the shocked-ffuid region gets compressed. This is 
determined using the conservation of rest mass of the shocked and unshocked fluids, which 
is given by 



Pfs(rs)^fs(rs) — Po(i)^o(i)- (12) 

Equation [12] provides the width of the shocked-fluid regions in their respective comoving 
frames and is given by 



A fs(r S ) = A o( i) • (13) 

The time of crossing the forward or reverse region by the FS or RS, in the shocked-fluid 
frame, can then be obtained by 

4,fs(rs) = n, { ] ■ (I 4 ) 
C ^fs(rs) 



2.2. Emission Region Parameters 

The calculation of the emitted radiation is dependent on the values of the magnetic 
field (B'), minimum (7mm) an d maximum (7 max ) electron Lorentz factors (ELF), and the 
normalization factor Q'$ m (cm^V 1 ) of the electron injection function Q' e m -'(7, t)(cm~ 3 s _1 ). 
The power-law distribution of the injected relativistic electrons is given by 



QeMrs)^') = Qo,fi(rs)7' q f° r 7 m in,f s (rs) — l' — 7 max ,fs(rs) (^) 

where q' is the electron injection power-law index. 

In order to proceed with the calculation of emission, we first derive the values of these 
parameters from shock dynamics (see below) for each of the shocked regions. We calculate 
the radiative energy loss rates of electrons and corresponding photon emissivities for both 
regions, separately, and add their respective contributions to the observed spectra. 
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Some fraction of the bulk kinetic energy density of the shocked fluid is converted into 
the magnetic and electron energy density by the resulting shocks. We define a magnetic field 
partition parameter e' B such that 



TV 

u B,fs(rs) 
^fs(rs) 



(16) 



where £7f s (rs) is the kinetic energy density of the shocked fluid. This is obtained by using 
Eqns. [7] & and is given by 



^ife(rs) 



Po(i)° 2 (Tfsirs) — l) 



7 r fs(rs) + 1 

7-1 



(17) 



The quantity U# = B^, is J8tt is the magnetic energy density. Equation [TBI yields the 
value of the magnetic field, which is assumed to be randomly oriented in space and tangled 
in the jet plasma. The maximum ELF of the injected electrons is obtained by balancing the 
power gained from the accelerat i on, i.e. from gyro-resona nt processes, with the synchrotron 
losses (jdeJager fc Hardingill992l ; IChiang fc Dermerlll999l ). This condition yields 



Tmax,fs(rs) 



4.6 x 10' 



a' 



£> fs(rs) 



where a' < 1 is the electron acceleration rate parameter, defined through a' = t' gyi /t' acc . 
The quantities t' gyI and t' acc are, respectively, the gyration and acceleration timescales of an 
electron. The quantity B' {s ,, is in units of Gauss. 

The minimum ELF is obtained by using the equation 



^e,fs(rs) — £ e^fs(rs) ■ (19) 

Here > ^e,fs(rs) = <,fs(rs)(7L(rs)) m eC 2 is the electron energy density, < fs(rs) = C^f is the 
number density of non-thermal electrons injected by shocks into their respective emission 
regions, (' e is the fraction of electrons that is accelerated behind the shock fronts into the 
power-law distribution given by Eqn. [TSJ and e' e is the electron energy partition parameter. 
The average electron energy ( 7 ^ (rs) ) for the energy range 7^ n ,f s(rs ) < l' is (rs) < 7L x , fs(rs ) is 
given by 
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'*^max,fs(rs) 

/ n e,fs(rs) (7f s (rs))7fs(rs)^7f s ( rs ) 

y 

/ / \ 'min,fs(rs) 

<7fe(rs)> = — ~, • ( 2 °) 

"max,fs(rs) 

/ ra e,fs(rs)(7f s (rs))^7f s ( rs ) 

^min,fs(rs) 

Here «e,fs(r S )(7fe(rs)) = n o7f s ~(rs) is tne number density per unit l[ s{rs) and n' is the initial 
electron number density. As most of the internal kinetic energy of the shocked fluid is 
stored in the baryons, implying that e' c « 1, we can write U' is ,-> = n' p fs ( rs )(rf s(rs ) — l)m p c 2 . 
Substituting all of the above expressions along with equation [18] into equation \T9\ we obtain 
the following expressions for j' m - n f ,y 



Chl(7min,fs(rs)) - 

7min,fs(rs) 7max,fs(rs) 
m (7min i fs(rs)) + ~j l n (7max,fs(rs)) 








7 



min,fs(rs) 



7 



if q' 
if q' 



1 

2 



max,fs(rs) 



y (W) rv , s 

'min,fs(rs)V 'min,fs(rs) 



'max,fs(rs)V fmax,fs(rs) 



)(7 m 



- C q ) = if q' y£ 1 or 2 



(21) 



where C = 1837^(rj s(rs) - 1), and C q = 1837jef f(T' is{rs) - 1). Equation EJJ is solved 
numerically to obtain the value of 7 / liT 1 fs ^ rs ^. In the case of ^> 7^ in and q' > 2, the above 
expressions reduce to 

f q' — 2 s' 

7min,fs(rs) = ^37 — -^(r fs ( rs ) — 1) . (22) 



The value of Qofe(rs) (^') * s obtained by assuming that a fraction of the kinetic en- 
ergy stored in the shocke d plasma is transferred, per unit time, to non-thermal electrons 
(IBottcher &: Dermerll2010l ). This yields 



f' TV 

& e u fs(rs) 



x,fs(rs) 



cr,fs(rs) 



( 5e,fs(rs)(7fs(rs)) " t ')7fs(rs) m e c2 d7f s (rs) ■ 



,fs(rs) 



(23) 



Using Eqns. [23] and [15], we obtain the expression for (£') as 



( 5o,fs(rs)(^) 



£ c^fs(rs) 



(2-9') 



c 2 t / i '(2-q') 

e cr.fs(rs) ^ 'max,fs(rs) 'min,fs(rs) 



V (2 - q ? A 

'min t r ^ \ ' 



iiq'^2 



m c c 2 t' 



c fs(rs) 



cr,fs(rs) ^ 'max,fs(rs) ' 'min,fs(rs) 



H 

1 'Tt 



if q' 



(24) 
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3. Photon Escape Timescale 

For a cylindrical geometry, the volume and angle-averaged photon escape timescale is 
not the same as that of a spherical one. The average photon escape timescale is used in 
the numerical computation of the evolution of the particle and photon population inside the 
emission region. Here, we derive a semi-analytical expression of the escape timescale and 
assume the cylindrical emission region to be Compton thin. 

We consider three possible directions by which a photon can escape from a cylindrical 
region: forward, sideways and backward. All quantities in this section and thereafter refer to 
the comoving frame only, unless otherwise indicated, hence the prime notation is no longer 
used. 

Figure [2] depicts the angle considerations for the three possible directions of escape for 
a photon from a cylindrical region of height (width) h and radius R. As stated earlier, the 
z-axis is along the axis of the cylinder. The volume- averaged photon escape timescale can 
then be written as 



R h 2tt 

(£ P h,esc) v = I I I (Wsc {r, z))rd(f)dzdr , (25) 
II ii ii 




V C yl 

where the angle-averaged escape timescale is written as 



27T +1 



(Wsc(r, z)) = ^- J J Z csc (/i, $; r, z)dfi d$ (26) 

o -1 

Here, l esc is the distance required to escape from the region in any direction, and fi = cos 9. 
We have solved the above integral semi-analytically to obtain the final expression for the 
volume and angle-averaged photon escape timescale for a cylindrical region, which is given in 
Eqn. [271 The intermediate steps needed to obtain the final expression are given in Appendix 
IB1 of this paper. 



2ir 1 1 



7 7 7 p n p 




2tt 1 

R 




, . xk arcsin = ) dxd<t> — 

vrc J J \y/k 2 + o 
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Fig. 2. — Illustration of the three possible directions of escape for a photon from a cylindrical 
region. The quantities l esc+ and l esc -, respectively, refer to the escape path lengths for a 
photon in forward and backward directions, and 8 + and tt — 6_ are the corresponding angles. 



2tt 1 



R 2 



II 



xk 2 In ( 1 + — ) dx d(j) , 




(27) 



2nhc 



o o 



where a = h/R and k = yl — (xsin0) 2 



— x cos 0. 



Using the geometry of a cylinder, the volume-averaged probability of escape for a photon 
in the forward direction is Pf w( j = 7iR 2 /(2nR(R + h)). The volume- averaged probability for 
a photon to go in the backward direction is the same as its counterpart in the forward 
direction, hence Pback = -Pfwd- Similarly, the volume-averaged probability of escape in the 
sideways direction is P s ;de = 2irRh/(2TcR(R + h)) = h/(R + h). These probabilities have 
been used to calculate the appropriate volume-averaged photon escape rates in the three 
directions, as described in the next section. 



The emission region parameters and the average photon escape timescale form the neces- 
sary set of quantities required to calculate radiative energy loss rates and photon emissivities. 
They are needed to follow the evolution of electron and photon populations inside the emis- 
sion region in a self-consistent and time-dependent manner along the length of the region. 
Since the same set of equations is applicable to both the forward and reverse emission regions, 
we do not distinguish between the various quantities used in the equations. 

We consider instantaneous acceleration of relativistic particles behind the shock front. 
The accelerated particles are injected into the region with a single-power-law distribution 
described by Eqn. [TSJ The acceleration of particles and the longitudinal expansion of the 
emission region continues until the shocks reach their respective boundaries of the merged 
shell. As the shocks propagate in their region, energetic electrons continue to be produced 
at the shock front. Owing to energy losses, only less energetic electrons will be found at 
progressively larger distances from the shock fronts. This creates an energy gradient within 
the emission region. 

Once the acceleration sets in, particles begin to lose their energy via synchrotron and 
SSC processes. The time- dependent evolution of the electron and photon populations, in 
each of the emission regions, is followed by using the equations 



4. 



Numerical Method 



dt 



97 [[it 



) 



loss 



n e ("/,t) +Q e (j,t) 




(28) 
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and 



dn ph (e,t) 



^ph,em\ 



Tlph,abs 



n ph (e,t) 



(29) 



^ph,esc 

Here, {d'j / dt)\ oss is the electron energy loss rate, due to synchrotron and SSC emission, 
Q e {l-,t) is the sum of external injection and intrinsic 7 — 7 pair production rate, and 
t e ,esc = r/R/c is the electron escape time scale with r\ being the escape parameter. The 
quantities n p h,em(e,£) and n p h,abs(e, t) are the photon emission and absorption rates cor- 
responding, respectively, to the electrons' radiative losses and gains, e = hv/m e c 2 is the 
dimensionless photon energy, and t p h,esc is the volume and angle-averaged photon escape 
timescale for a cylindrical emission region (Eqn. |27|) . In order to obtain the temporal 
evolution of the electron distribution in an emission region, equation is discretized and 
rearranged in such a way that it forms a tridiago nal matrix, which is then solved numerically 
Jchiaberge fc Ghisellinilll999l : IPress et al.lll992h . 



Our current approach simulates the early phase of 7-ray production. During this phase, 
the radiative cooling dominates over the adiabatic cooling and the emission region remains 
optically thick out to < GHz radio frequencies. As a result, the simulated radio flux is well 
below the observed radio data. We do not model the phase in which the jet components 
gradually become transparent to radio frequencies because that would introduce several 
additional and poorly constrained parameters. 



The synchrotron loss rate is calculated from (jRybicki fc Lightmarull979l ) 
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syn 
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5 2 7 2 



57rm P c z 



(30) 



where <tt = 6.65 x 10 -25 cm 2 is the Thomson cross-section, and B is obtained from Eqn. [16j 
The synchrotron photon prod uction rate per unit volume in the energy interval [e, e + de] is 
calculated using the formula (ICrusius &: Schlickeiserlll986[ ) 
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where x 



Airm e cu 
3eB-y 2 



A numerically simplified version of the synchrotron emissivity (described 



in Appendix [C]) has been used in our simulations to save CPU time. 

We have incorporated the effects of SSA on the sync hrotron spectrum in our ca lculation. 
The SSA optical depth is calculated using the formula ( jRybicki fc Lightmanl 119791 ) 
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(32) 



where / p h,esc = £ P h,esc c is the mean path length traversed by a photon escaping from its point 
of origin inside a cylindrical region. The synchrotron emission from a cylindrical region is 
then obtained using the expression 



exp(-rssA)) 

r SSA 



7lR 2 h . 



(33) 



The electron energy loss rate due to scattering an isotropic, monochromatic radiation 
field of photon energy e and density n p h (cm -3 ), including the effects of the Klein-Nishina 
scatte ring cross section, is calculated using equation 32 of IBottcher. Mause &: Schlickeiser 
(119970 . The pho ton product i on rate per unit volume of isotropic S SC emission in the region 
is calculated by (I Jones! 1 1968c IBottcher. Mause fc Schlickeiser! 119971 ) 



ri S sc(e s ,fis) = j djn e (-f) J den ph (e)g(e s ,e,-f) , (34) 
1 o 

where e s is the energy of the scattered photon, and n ph (e) is the radiation field at photon 
energy e available for SSC scattering in the region. The quantity q(e Hl e,j) is a function that 
describes the probability of Compton scattering, and is given by ( jJoneslll968! ) 



3ccr T ( 47 2 e s _ i 



^(e s ,e,7) 



167 4 e 



3cctx 
47 2 e 



for ^2 < e s < e 



(35) 



2gln(g) + (l + 2g)(l-g) + 



(4e 7 g) 2 (l-g) 



for e < e s < 



4e 7 2 
l+4e 7 ' 



where q = e s /4ej 2 (l —e s /j). The radiation field available for SSC scattering includes all 
contributions from the photon field of the previous time step. This accounts for SSC scat- 
tering to arbitrarily higher orders. The optical de pth for a 7-ray photon of energy ei du e to 
7 — 7 absorption is calculated using the formula (iBottcher. Mause fc Schlickeiser! 119971 ) 



+ 1 00 

27r/ ph;CSC J dn{l-y) J deo- 77 (ei,e,yu)ri ph (e, Q) , 



(36) 
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where n p h is the photon density present in the emission region that provides the photon 
field for 7 — 7 absorption, and <r 77 is the pair production cross-section. The subsequent 
pair production rate is obtain ed using the exact and analytic solution, given in Eqn. (26) of 



Bottcher fc Schlickeiserl (119971 ). The contribution from pair production is added to Eqn. 
at every time step of the simulation. Using the radiative transfer equation, the high-energy 
emission that is able to escape from the region can be calculated as 

Nifci^s, a) = n ssc (e s , a) (1 ~ eXP( ~ T77) W fr . (37) 

T 77 

The temporal evolution of the photon population inside the emission region, for synchrotron 
and SSC emission, is followed using equation [29j The time step used in the simulations is a 
fraction of the minimum time scale among the relevant ones (cooling, electron and photon 
escape, injection, and shock crossing for as long as the shocks are in the region) from both 
regions. The simulation time step is taken to be common for both emission regions. 



4.1. The Slice Scheme 



In order to reproduce the observed spectral variability patterns of a blazar, it is im- 
portant to consider the inhomogeneity in the particle as well as photon density throughout 
the emission region. Observations in the optical and higher energy bands indicate that the 
acceler ation and /or cooling timescal es can be shorter than the light crossing time for the 
region ( Chiaberge Sz Ghisellini 19991 ). The highly energetic electrons that have been accel- 
erated freshly at the shock front evolve on a cooling timescale shorter than the light crossing 
time, farther from the shock front those ultra-high energy electrons have had time to cool 
significantly. This induces an energy gradient in the electron population within the region 
and creates an inhomogeneity in the particle as well as photon energy density throughout 
the emission region. As a result, the observer se es a combination of yarious spectra b eing 
produced in different parts of the region (see also lSokolov et al.l 12004 ; iGraff et al.l 120081 ) . 



In our model, we have incorporated the inhomogeneity in the particle and photon energy 
densities by dividing both emission regions into multiple slices, each of width h z and radius 
R. As shown in Figure |3l for simplicity we start with a jet that is devoid of relativistic 
particles. The first population of ultrarelativistic particles is injected by the shocks very 
close to the CD. 

The slice closest to the CD in the forward (reverse) emission region will have the forward 
(reverse) shock propagating through it first. As the shock advances, it injects particles 
into the slice with a power-law distribution that is dictated by the shock parameters. The 
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Fig. 3. — Illustration of the slice scheme used to calculate the integrated radiation spectrum 
resulting from the emission regions shown in Fig. [TJ The entire forward and reverse emission 
regions are divided into iVf s and iV rs slices, each of width h z . The quantities c x n p h,f„d, c x 
%>h,back, andc x n p h )S ide are the photon fluxes (photons per unit time per unit area of the 
slice) in the forward, backward, and sideways directions delineating the transfer of photons 
in between the slices throughout the forward and reverse emission regions. 
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normalization factor, Q™ 3 is calculated for every slice using equation [2H and replacing C/fe( rs ) 

by t4(rs),zn = 7^7 and £cr,f s (rs) by t cr , fs ( rs ), zn = ^t 2 — ■ The quantities described in g] are 
calculated for each slice. As the shock enters a new slice, it gradually populates that slice 
and the entire process of injection, acceleration and cooling is repeated. For sake of brevity, 
we refer to this phase of a slice as the acceleration phase in the following text. The injection 
and acceleration stops in the previous slice, which the shock has just left, and only the 
cooling continues. We refer to this phase of the slice as the cooling phase. 

The photon density of a slice, along with Pf w( j, -Pback, and P si d e (see $3]), are used to 
calculate photon escape rates in all three directions, according to 

dn phiiwd (e,tt) _ n ph (e,Q) 

dt ^ph,esc 

and 



fwd 



(38) 



dn ph:Sidc (e,tt) _ n ph (e,tt) 

t, — , -T side j I"" J 

w£ £ph,esc 

with driph,back(e, ^) = ^%>h,fwd(e, The sideways propagating photon density becomes part 
of the observed spectrum originating from that slice. The forward and backward propagating 
photon densities are used to provide feedback to the neighbouring slices. The density of 
photons escaping from the region is subtracted from the total density of that slice for the 
current time step and the remainder is used for calculations in the next time step. We point 
out here that even if the volume and angle-averaged photon escape time scale for a slice 
turns out to be greater than the simulation time step, on average one would still expect 
some photons to escape out of the slice, as soon as they are created, if they happen to 
be produced close to the boundary. In that respect, our assumption that some percentage 
of photons do escape out of the region is justified and can be used to provide feedback to 
adjacent slices without having to track their movement individually at every time step. 

The process of providing photon feedback from a particular slice to its adjacent slices, 
throughout the emission region (forward and reverse), at every time step accounts for the 
inhomogeneity in the photon density in the jet. In order to accurately calculate the SSC 
process, the non-local retarded-time SSC losses due to photons produced in other parts of 
the source need to be taken into account. These are very important for SS C dominated 



sources, such as TeV blazars, and thus cannot be ignored (IGraff et al.l 120081 ) . Our multi 



slice scheme with radiation feedback, lets us inherently addr ess this issue without having to 



calcu l ate the photon density rat es at retarded times (see also lSokolov et al.ll2004j ; iGraff et al. 



20081 ; iBottcher fc Dermerl 120101 ). along the length of the emission region. We note that our 
approach does not separately address the retarded-time photon field in the radial direction 
as calculating the SSC process in 2D is computationlly more intensive and beyond the scope 
of this work. 
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5. Light Travel Time 



Since photons travel with a finite speed the time delays should be incorporated in the 
calculation of arrival time of the radiation in the observer's frame. Different distances are 
covered by photons originating from different parts of the source but reaching the observer 
at the same time. In the case of blazars, the axis of the jet makes an angle 9* obs with the 
observer's line of sight. Since the emission region moves towards us, we need to take into 
account only the radiation that is emitted in the direction of 8* hs . Hence to calculate the time 
delay of the emitted photons, we consider the sideways direction for both emission regions, 
and the forward direction for only the forward emission region. In this section, starred 
quantities refer to the observer's frame, whereas primed quantities indicate the comoving 
frame of the emission region. 

In case of the sideways direction, the time delay equation for the i th slice, of width h' z , 
in the emission region is given by 



A *sidc,a(*a) = • (40) 

C 

Here N tot = Nf s + N IB is the total number of slices in the entire emission region, and 9' is the 
jet viewing angle given by 

cose'= COs9 ^~^ . (41) 
1- fact* 6^ 1 ; 

The subscript 'a' stands for the fs (forward) or rs (reverse) emission region. The value of 
iVb = when the time delay for the forward region is calculated and iVb = Nf B when the 
same is calculated for the reverse region. 

The distance calculation using Eqn SD] becomes effective only after the shocks leave a 
given slice. Until then, the precise location of the shock fronts is used to calculate the time 
delay in the sideways direction: 

4 = N h ti z ,{ 9 - cf3 fs t' (42) 

and 

4 = Nf B h z ,is + cfrst' , (43) 

where t' is the total simulation time elapsed at the current time step. The above two 
equations incorporate the fact that the space in front of the shocks is still devoid of ultra- 
relativistic particles and hence is inactive. In the forward direction, the time delay from the 
forward emission region is calculated as 
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, , IL cos 6' 

AC P ,fs = JL — • (44) 



We carry out the forward time-delay calculation for as long as the shock is inside the 
forward emission region. Once the shock leaves the region, only the forward radiation origi- 
nating from the slice closest to the observer's line of sight is observed. The forward radiation 
from other slices provides feedback to adjacent slices. 

All time delays are calculated with respect to the closest slice to the observer. In 
addition to the above mentioned time delays, the beaming of the radiation as observed 
in the observer's frame also needs to be included. This effect is incorporated using the 
Doppler boosting factor, D = [r s h(l — (3 s h cos 6** bs )] _1 , that connects the comoving frame of 
the emission region to the observer's frame. The total time-delay in the reception of photons 
from either direction, in the observer's frame, for the i th slice of an emission region is given 
by 

1 + Z 

Cid°c,fs(rs);up,fs = + ^sidc,fs(rs); uPj fe(i)) , (45) 

where Z is the redshift of the source. 

The frequencies have been transformed to the observer's frame by multiplying by D, 
whereas the energy fluxes per unit frequency have been transformed by multiplying by D 3 
such that (vF u )* in Jy Hz is given by 

(vF u )* = ■ (46) 

Here N' h is the sum of the comoving photon escape rates in the forward and/or sideways 
direction, from both emission regions, at the current time step, and is the luminosity 
distance to the source. 

In our model, we can accomodate the case of cos#* bs < /3 s h- This would imply that 
in the comoving frame of the emission region, the observer is located outside the cone of 
superluminal motion and is looking at the region from behind. As a result, the backside of 
the region would be visible to the observer sooner than the front side. Thus, in this case, the 
observer would see the sideways and backward emission of the reverse region before seeing 
the sideways emission from the forward emission region. Since the emission in the direction of 
9* hs needs to be considered, no radiation in the forward direction from the forward emission 
region would be visible to the observer. 



-22 - 



6. Parameter Study 

The combination of various input parameters dictates the overall evolution of the radi- 
ation spectrum and lightcurves in a simulation. This makes it important to understand the 
effects of varying such parameters in order to reproduce the observed features of a source. 
These parameters can be broadly categorized into shock, emission region, and jet parameters. 
In what follows, the unprimed quantities refer to the AGN frame, primed quantities to the 
comoving frame, and starred quantities indicate the observer's frame. The shock parameters 
that affect the overall evolution of the SED are the kinetic luminosity, L w , the total duration 
of the ejection event, i w , the mass of the outer shell, M Q , and the widths and BLFs of the 
inner and outer shells, Aj & T; and A Q & T . The emission region parameters that have an 
impact on the evolution of the radiation spectrum are the equipartition parameters, e' e & 
e' B (see §2.2p . the particle acceleration fraction, Q, the electron acceleration rate parameter, 
a', and the particle injection index, q' . Finally, the jet parameters that play a role in the 
evolution of the spectrum are the radius of the slices and the jet, R', and the orientation 
angle of the jet, 6>* bs . 

We have carried out approximately 25 simulations to test and study the effects of vary- 
ing the values of each of these parameters on the time-averaged simulated SEDs and their 
corresponding lightcurves. For all our simulations, the flux values are calculated for the fre- 
quency range v' = (10 8 — 10 26 ) Hz and electron energy distribution (EED) range 7' = 10 — 10 8 
with both ranges divided into 150 grid points. The entire emission region is divided into 
100 slices with 50 slices in the forward and 50 in the reverse shock region. We constrain the 
values of some of the parameters analytically to make sure that unphysical values are not 
used in the simulations. The acceleration time scale of the highest energy electron should be 
less than their corresponding synchrotron cooling time scale. The following condition is used 
to place an upper limit on the value of 7 max for particles in the both forward and reverse 
shock emission regions: 

^ £ S • (47) 

where e = 4.8 x 10~ 10 esu is the electron's charge in cgs units. The value of 7^. should not 
exceed the maximum value of the EED range, mentioned above. In case of ~ ',;„■ the value 



should not fall below the lowest value of the EED range. As pointed out by iMimica et al. 



( 120041 ) . the Larmor radius, r' L , of the fastest moving (highest-energy) electron should be 
smaller than the slice width so that the magnetic field strength exceeds a certain minimum 
value and the shock acceleration can take place. This condition is used to make sure that 
the number of slices for both regions is selected in such a way that 

< 1 . (48) 



r'r m 



hi eB' hi 



-23 - 



In order to carry out the parameter study, the value of each of the shock, emission region, 
and jet parameters is varied twice. We study the effects of variation in terms of changes in 
the flux level, spectral hardness (SH) of the simulated SED, and simulated lightcurves for 
different photon energies. Table [TJ shows the values of the base set (run 1) parameters used 
to obtain the baseline model. 

Table [2] shows the values of each of the parameters that are varied in the rest of the 
simulations. We study the effect on the simulated SED and lightcurves with respect to that 
of the baseline model. 

Simulations 2 and 3 involve changing the value of M Q as well in order to satisfy Eqn. 
[TJ and obtain a non-neg ative value of M h where M; = Lwtw r {A i° r ° c2) . On the other hand, 
changing the value of M Q only could be carried out just once because changing it in the 
other direction resulted in Mj > M Q , while we have assumed that the inner shell is the faster 
moving one and thus should have a relatively lower mass than the slower moving outer shell. 

7. Results 

Figure H] shows the resultant time-integrated SED and lightcurves of the baseline model 
for a generic blazar source over an 8-day period. The input parameters chosen for the base 
set result in a value for the BLF of the emission region as T s h = 14.9, with the values of BLF 
of the forward and reverse shocks being T' {s = 1.08 and T' rs = 1.14. A magnetic field value 
of B' = 2.51 G is obtained for both the forward and reverse emission regions. The values 
of 7^ in and 7^ ax resulting from the input parameters, and evaluated using equations [5T] & 
[TBI are obtained as 7mi n ,f s = 2.18 x 10 3 , 7^ in rs = 3.74 x 10 3 , and 7max ; fs,rs = 8-31 x 10 4 . The 
derived total width of the forward and reverse emission region is A' {s = 8.19 x 10 15 cm and 
A' TS = 9.93 x 10 15 cm, with the shock crossing time for each of the emission regions being 
t' cr is = 7.20 x 10 5 s and t' CTTS = 6.94 x 10 5 s. 

The time-integrated SED of the base set shows that the synchrotron flux peaks in the 
optical-near UV at z/* yn = 8.59 x 10 14 Hz, whereas the high -energy SSC component peaks in 
the MeV regime at ^gg^ — 8.72 x 10 21 Hz. The synchrotron to SSC transition takes place 
in the hard X-rays at z/ t * urn = 2.96 x 10 17 Hz, implying that synchrotron photons extend into 
the soft X-ray regime for this generic blazar source. The derived Compton dominance factor 
(CDF), defined as CDF = uF* ssc.pcak/^* syn.peak^ ig 2 .go. The photon spectral index in 
the X-ray (2 - 10 keV) range is found to be a^-iokev = 0.549, whereas in the Fermi range 
(calculated at 10 GeV) it is a* WGeW = 2.08. 

The right side of Figure H] shows the simulated lightcurves of the baseline model for the 
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Table 1. Parameter list of run 1 used to obtain the baseline model. 



Parameter Symbol Value 



Kinetic Luminosity 


L w 


10 47 crg/s 


Event Duration 


tw 


10 7 s 


Outer Shell Mass 


Mo 


4 X 10 31 g 


Inner Shell BLF 


r t 


25 


Outer Shell BLF 


r D 


10 


Inner Shell Width 


A t 


3 X 10 15 cm 


Outer Shell Width 




6 X 10 15 cm 


Electron Energy Equipartition Parameter 


< 


5 x lO" 1 


Magnetic Energy Equipartition Parameter 




2 x 10~ 3 


Fraction of Accelerated Electrons 


a 


2 x 10~ 2 


Acceleration Timescale Parameter 


a' 


8 x 10- 6 


Particle Injection Index 


q' 


3.4 


Slice/ Jet Radius 


K 


3 x 10 16 cm 


Observer Frame Observing Angle 


%bs 


3.15 deg 


Redshift 


Z* 


0.306 



Table 2. Parameter list for other simulations. 



Run # 




Parameter Value 


2 




L w = 5 X 10 47 ergs/s 


3 




L w = 1 X 10 46 ergs/s 


4 




M = 1 x 10 32 g 


5 




Ti = 18 


6 




Ti = 30 


7 

8 




r D = 5 
r D = 14 


9 


A, 


= 3 x 10 16 & A = 6 x 10 16 cm 


10 


A, 


= 6 x 10 14 & A D = 9 x 10 14 cm 


11 




4 = 9.5 x lO" 1 


12 




4 = 1 x 10~ 2 


13 




s' B = lx 10-' 2 


14 




e' B = 2x 10-'' 


15 




C e = 9.5 x 10" 1 


16 




C' e = 2 x 10- 3 


17 




a' = 6 x 10~ 5 


18 




a' = 8 x 10~ 8 


19 




q' = 4.6 


20 




q' = 2.0 


21 




R' z =3x 10 17 cm 


22 




R' z = 1.4 X 10 16 cm 


23 




e :b s =3-8dcg 


24 




^ bs = 1 - 1 5deg 


25 




6»* bs = 10.0 deg 
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energy bands in the optical (R), X-ray (10 keV), high-energy (1 MeV), and VHE regime 
(1 TeV). The flux value for each energy band in all the graphs has been normalized with 
respect to the maximum flux value of that energy band. 

The lightcurves indicate that the optical flux continues to rise as long as the shocks are 
present in the emission region. It peaks at t* cak = 54 ks. Once the shocks break out of the 
system, cooling dominates and the pulse begins to decline steadily. Thus, the duration for 
which the optical pulse lasts at its half maximum (FWHM) is 0.91 days in the observer's 
frame. Optical synchrotron photons are produced from higher energy electrons, which cool on 
a time scale shorter than the dynamical time scale of photons {t' dyn = l' csc /c) for a particular 
slice. Here l' esc is the comoving average escape length for a photon originating anywhere in a 
slice. As a result, the distribution of such electrons remains steady during the acceleration 
phase of a slice. This implies that the total optical flux includes a contribution from the 
accelerating slices and the observer sees an increase in the flux due to an increase in the 
volume of the emission region. Once the shocks exit the respective emission regions and the 
cooling phase dominates the region, the optical flux decreases immediately. This makes the 
rising and decaying phases of the pulse almost equal leading to a quasi-symmetrical profile 
of the optical pulse. 

On the other hand, the X-ray flux continues to rise even after the shocks exit the 
emission region. This happens because of two reasons: the 10 keV photons are produced 
from upscattering of lower energy (near IR) synchrotron photons off lower energy electrons. 
Such electrons remain in the system for a longer time due to slow cooling. As a result, the 
rise in the corresponding intensity is slow. Also, since no more electrons are added into 
the system during the cooling phase there is a continued increase of late-arriving photons 
at the sites of scattering due to the internal light crossing time effect. As a result, initially 
we receive the emission from the distribution that starts to cool in the back slices while the 
front slices are still accelerating. Once the cooling phase takes over the emission region, an 
important contribution comes from the front slices as the back slices have a higher number 
of cooled off lowest energy electrons. As a result, for a SSC dominated radiation of v ~ 10 16 
Hz, the pulse experiences delayed peaking and decays with a long and slow tail compared to 
its rise. This explains why the 10 keV pulse has an asymmetrical profile and peaks last at 
t* eak = 81 ks and persists for the longest time compared to the pulses at other frequencies, 
with a FWHM of 1.8 days (in the observer's frame). Dividing the emission region into various 
slices increases the spatial resolution of the region. As a result of this enhanced resolution, 
when the shocks exit the system and pulses of higher synchrotron and SSC frequencies die 
due to the reason explained above, the 10 keV pulse, which builds up slowly over time bears 
the signature of this exit. The kink seen in the lightcurve of the 10 keV photons is the 
signature of shocks completely breaking out of the system with the forward shock - the last 
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shock to come out of the region - exiting the forward emission region at t* r = 5.27 x 10 4 s. 

The 1 MeV photons are a result of Compton upscattering of higher energy optical 
photons off low energy electrons, making the pulse peak at the same time as R-band photons 
at t* k = 54 ks. At all energies, the pulse rises sharply and declines comparatively gradually. 
This implies that, for the reasons explained in the case of the R-band pulse profile, the profile 
of the 1 MeV pulse is also quasi-symmetrical. 

The very high energy (VHE) photons (1 TeV) result from the SSC scattering of higher 
energy (soft X-rays) synchrotron photons off higher energy electrons. As these electrons 
cool at a very fast rate compared to the electrons emitting at optical frequencies, the VHE 
photons are produced only during the acceleration phase. This is why the corresponding 
pulse has a sharper rise and peaks sooner than its lower energy counterparts, at t* eak = 45 
ks. The pulse undergoes a sharp drop in the 1 TeV energy regime exactly at the time 
the shocks completely leave the system, thereby having the shortest pulse duration with a 
FWHM of 0.81 days and a nearly symmetrical profile. 

For all simulations, including the baseline model, the average SED has been evaluated 
over an integration time of 9000 s. This facilitates comparison with observations, which 
usually require extended exposure times. Table |3] lists the effects of changing the values 
of the input parameters (according to Table |2]) on the flux, the peak synchrotron and SSC 
frequencies (z/* yn & ^ssc' respectively), the trough frequency (^ t * urn ) indicating the transition 
from the synchrotron to SSC component, the SH represented by energy spectral indices a* 
in the 2-10 keV range and at 10 GeV, and the CDF of the baseline model's time-integrated 
SED. 

Table H] parametrizes the simulated lightcurves in terms of the time of the last shock to 
break out of the region, t* T , the time at which the pulse at a given energy peaks, t* eak , and 
the FWHM (indicated as W* in the table) of the pulse at a given energy. 

The simulated time-integrated SEDs and the corresponding lightcurves are affected in 
several ways as input parameters are varied. The variation directly changes the values of 
some of the physical quantities like r' fs , T' rs , relative velocity of the shells, internal energy 
of the shocks, magnetic field, 7^ in & 7max> particle density, and size of the region, which 
affect the overall flux, SH, and CDF of the spectrum. This also directly alters the time at 
which a pulse at a given energy peaks, as well as its FWHM. We describe below some of the 
dominant effects observed on the simulated SED and lightcurves according to the changes 
in the values of the above mentioned physical quantities. 

A) An increase in the internal energy of the shocks boosts the flux level of the spectrum. 
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Table 3. Parameter study of shock, emission-region, and jet parameters on the 

time-integrated SED. 



Input Parameter 


Value 


vF* 


"syn [ Hz l 


"tain l Hz l 


^ssc 


/ ( T'l 771* 


C IT* /O 1fl 1 t\ 
ori [Z-W keV ) 


O TJ* /in (^„1A 

brl (1U Oev ) 


Baseline 






8.59 x 10 14 


2.96 x 10 17 i 


3.72 xlO 21 


2.90 


rvn mi™v = 0.549 

^ Z — 1UKC V J J J 


OiTnn^r — 2.08 




t 


t 


t 


t 

1 


t 


3.09 


I (\n , nl -tr — 6S2 


f rv-i n p„ v = 2 04 




4- 


4- 


4 


t 


4- 


2.20 


t "2-iOkcV = 0.444 


4 "lOGcV = 2.09 


M a 


t 


4- 


t 


t 


t 


4.84 


4- "2-iOkcV = 0.847 


t aioGcV = 1-75 


r< 


t 


t 


t 


t 


t 


2.86 


4- "2-iOkcV = 0.824 


t aioGoV = 1-92 




4- 


; 


4- 


4- 


4- 


2.65 


t "2-lOkcV = 0.441 


4 "lOGcV = 2.23 


r D 


t 


4- 


4- 


4- 


4- 


2.34 


t "2-iOkcV = 0.377 
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Note. — The slope of the low-energy component of the time-integrated SED obtained from run 1 in the X-ray range 
(2-10 keV) is c*2— iokeV — 0.549 and that of the high-energy component in the Fermi range (10 GeV) is ^inrjev — 2.08. 
The time-integrated low-energy component of the baseline model peaks in the optical-UV range at ^g yn = 8.59 X 10 14 
Hz. The high-energy component peaks in the MeV range at fg gc = 8.72 X 10 21 with the transition from synchrotron to 
SSC component taking place in the X-rays at v* UIn = 2.96 X 10 17 . The f implies an increase in the value of a particular 
quantity, J. refers to a decrease in the value, and — stands for no change in the value. The abbreviation CDF stands 
for Compton Dominance Factor, defined in J7] and SH stands for Spectral Hardness, defined in 
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Table 4. Parameter study of shock, emission-region, and jet parameters on simulated 

lightcurves. 



R-Band 10 keV 1 MeV 1 TeV 



Input Parameter 


Value 


t* r [days] 




W [days] 


*pcak M 


W* [days] 


*peak H 


W [days] 


pcaK L J 


W* [days] 


Baseline 




0.61 


5 1 


0.91 


81 


1.8 


5 1 


1.1 


45 


0.81 


L w 


t 


0.61 


45 


0.60 


54 


1.1 


45 


0.78 


45 


0.56 




4 


0.61 


45 


1.1 


180 


3.2 


99 


2.3 


45 


0.56 


Mo 


t 


1.6 


63 


1.6 


162 


2.1 


144 


2.2 


54 


1.5 


r, 


t 


0.55 


45 


0.78 


72 


1.8 


54 


1.1 


45 


0.95 




4 


0.83 


72 


0.97 


108 


1.9 


72 


1.3 


63 


0.85 


r 


t 


1.3 


72 


1.5 


144 


2.3 


99 


1.7 


72 


1.2 




4 


0.37 


18 


0.34 


36 


1.6 


18 


0.68 


18 


1.4 




t 


6.1 


522 


5.5 


501 


4.4 


195 


5.0 


522 


5.5 




4 


0.11 


9 


0.27 


18 


0.97 


9 


0.40 





0.26 


4 


t 


0.61 


54 


0.85 


72 


1.8 


54 


1.1 


45 


0.76 




4 


0.61 


45 


0.61 


45 


2.2 


144 


3.5 


45 


0.50 




t 


0.61 


45 


0.56 


15 


1.3 


45 


0.80 


45 


0.54 




4- 
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63 


1.39 
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2.44 


99 


1.7 


54 


1.2 


c 


t 


0.61 
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3.3 


162 


3.4 
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2.4 


45 


0.88 




4 


0.61 


45 


0.64 


45 


0.54 


15 


1.7 


45 


0.54 


a' 


t 


0.61 


54 


0.91 


81 


1.8 


54 


1.1 


45 


0.81 




4 


0.61 


15 


0.62 


72 


1.9 


45 


0.94 


45 


0.58 


q' 


t 


0.61 


54 


0.91 


72 
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54 


1.0 


45 


0.88 




4 


0.61 


45 


0.88 


45 


1.1 


54 


1.3 


45 


0.63 




t 


0.61 


117 


7.4 


1332 


28 


540 


18 


54 


0.60 




4 


0.61 


45 


0.55 


45 


0.68 


45 


0.59 


45 


0.55 


n* 

%bs 


t 


0.72 


54 


1.1 


90 


2.1 


63 


1.5 


54 


0.83 




4 


0.40 


36 


0.70 


54 


1.3 


36 


0.77 


27 


0.81 


outside beaming cone 


t 


2.8 


23 1 


1.9 


432 


5.7 


252 


5.1 


234 


2.0 



Note. — The lightcurves of the baseline model are parametrized to study the effects of variation of various input parameters, 
listed in Table [2] We parametrize these lightcurves according to the shock crossing time, t* T1 time of peaking of a pulse at a given 
energy, £* eak , and full width at half maximum (FWHM), indicated as W* , of a pulse at a particular energy. 
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B) An increase in the density of injected electrons increases the overall flux of the 
broadband spectrum. It also increases the photon number density in the region. This 
implies rapid radiative cooling for electrons due to synchrotron and SSC processes, and a 
higher value of CDF due to an increased level of SSC flux. In terms of lightcurve profiles, 
this translates into acceleration and cooling setting in sooner, thereby making the pulse at all 
energies peak sooner and last for a shorter time compared to the case where the acceleration 
and cooling set in more gradually. 

C) An increase in the value of relative velocity of the two shells, along with the value 
of their individual shock propagation Lorentz factors (rj- s & T' rs ), increases the acceleration 
efficiency of the collision. As a result, the overall flux of the spectrum rises. In terms of 
lightcurve profiles, due to efficient acceleration the particles are energized sooner and start 
to cool down faster. Thus, the higher the acceleration efficiency, the faster is the particle 
acceleration and the subsequent cooling. Hence, pulses at all energies peak sooner and last 
for a shorter duration of time, except for the VHE pulse. The FWHM of the VHE pulse lasts 
for a longer time, compared to that of the baseline model, because of efficient acceleration 
that energizes more particles to higher energies. As a result, high-energy electrons capable 
of producing TeV electrons are available for a longer time, and the FWHM of the TeV pulse 
increases. 

D) An increase in the value of 7^ nin implies more higher energy electrons being injected 
into the system. The resultant spectrum shifts toward higher frequencies. The shift in peak 
(z/* yn & ^ssc) an d trough (v* UTrL ) frequencies makes the spectrum softer in the X-ray (2 - 
10 keV) but harder in the Fermi range (10 GeV). Usually the X-ray range forms the lower 
energy part of the SSC spectrum while the Fermi range is around or beyond the SSC peak. 
However, if v* y increases, the high-frequency end of the synchrotron component may extend 
into the 2-10 keV regime. As a result, the corresponding value of c^-iokcv increases, 
indicating a softer X-ray range spectrum. A shift of z/g SC to higher frequencies will harden 
the Fermi spectrum. 

E) A much lower value of 7^ in implies the presence of very low energy electrons in the 
system. This shifts the spectrum leftwards, toward lower frequencies, and higher order SSC 
components become prominent. This happens because the Thomson depth of the region 
increases and Klein-Nishina effects no longer suppress the higher-order SSC components. 
In terms of lightcurve profiles, this implies that acceleration and hence cooling set in very 
gradually for low energy electrons. On the other hand, higher energy electrons do not last 
for a long time in the system. As a result, photons resulting directly from these electrons 
are produced very early on and last for a very short time. 

F) The magnetic field in the emission region is responsible for synchrotron radiation. An 
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increase in the value of the magnetic field decreases the synchrotron cooling time scale, which 
makes the synchrotron mechanism the primary mode of radiative cooling for the electrons. 
As a result, the synchrotron component of the spectrum becomes dominant while the photon 
flux due to SSC scattering decreases, thus reducing the resultant CDF value. The overall 
flux of the spectrum increases, but due to excessive radiative cooling the spectrum becomes 
softer, as it comes from a cooled population of electrons. In terms of lightcurve profiles, 
the shorter synchrotron cooling time scales lead to shorter pulses compared to the baseline 
model. 

G) A larger emission region implies a larger comoving volume, but lower values of the 
magnetic field, shocks' internal energies, and density of injected electrons in the region. As 
a result, the flux and CDF of the spectrum decrease following the explanation given in (A) 
and (B). In terms of lightcurve profiles, since the acceleration and radiative cooling set in 
gradually, the pulse at a given energy peaks later and has a higher FWHM relative to that 
of the baseline model. 

H) An increase in the value of 7m ax without any change in the magnetic field value implies 
a shorter synchrotron cooling timescale for the highest-energy electrons. This results in more 
synchrotron photons in the region, thereby boosting the flux level of the SSC emission. 

Subsections 17.11 - 17.31 briefly discuss the effects of varying various input parameters on 
the time-averaged simulated spectra and lightcurves of a generic blazar source. 

7.1. Shock Parameters 

Figures [5] - [7] depict the effects of varying shock parameters (L w , M a , T;, r o , and Aj 
& A ) on the simulated time-integrated SED (left side) and lightcurves (right side) of the 
baseline model (run 1). 

Increasing L w (run 2) increases the value of the internal energy of the shocks, the 
magnetic field, an d the density of the injected electrons in both emission regions. Due 
to the reasons cited in (A), (B), and (D) the overall flux of the spectrum rises (see Fig. EJ) 
and the SH decreases in the X-ray range but increases in the Fermi range (see Table |3]). 
The CDF value increases marginally due to the reason cited in (B). On the other hand, 
decreasing L w (run 3) from its original value results in a spectrally harder spectrum in the 
X-rays and a softer one in the Fermi range, along with an overall lower flux of the broadband 
time-integrated SED and a reduced CDF value. 

Making the outer shell heavier (\ M , run 4) makes the shell denser while making the 
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inner shell lighter and rarer. This combination of densities decreases the values of r s h and T' {s . 
The shocks' internal energies and the magnetic field values decrease along with the density 
of electrons injected into the forward emission region. This results in an overall decrease 
of the flux throughout the spectrum according to the reasons cited in (A) and (B). The 
7minfs vame decreases, whereas 7mm.rs increases due to an increase in r^. Consequently, a 
mixed population of lower and higher energy electrons arises in the emission region during the 
cooling phase, giving rise to a bumpy time- integrated SED. This also leads to a comparatively 
longer acceleration and cooling phase in the emission region. Since the X-ray range falls in 
the 3rd hump of the spectrum, which forms the lower energy part of the SSC component, the 
spectrum is softer (see Table [3]) in this case. On the other hand, an increase in 7^ in rs results 
in the shifting of u* yn , ^ t * urn , and z/g SC to higher frequencies. As a result, the spectrum is 
harder in the Fermi range, as explained in (D) and indicated by the value of ai GeV m Table 
[21 An increase in T' TS decreases the comoving width of the reverse emission region, thereby 
decreasing the comoving volume and increasing the comoving number density of electrons in 
the region. This leads to an increased contribution of SSC photons from the reverse emission 
region, thus slightly increasing the CDF value. 

The lightcurves of run 2 (L w f), 3 (L w 1), and 4 (M Q f; right side of Fig. E]) indicate a 
similar pattern as that of run 1. The same reasons used to explain the time of peaking of a 
pulse and its FWHM, at a given energy, for run 1 in JT]can be applied to runs 2, 3, and 4 to 
understand their pulse profiles. In case of run 2, pulses at all energies peak sooner compared 
to run 1, except for the TeV pulse. They also last for a shorter time due to the reasons cited 
in (B). The opposite holds true for run 3, although the optical pulse peaks sooner while the 
TeV pulse has a smaller FWHM (see Table HJ) compared to run 1. This happens because the 
R-band synchrotron photons and the TeV Compton upscattered photons are now a result 
of relatively higher energy electrons, in comparison to run 1. In the case of run 4, the pulse 
peaks later and lasts longer at all energies compared to that of run 1 (see Table 0J due to 
the reasons mentioned in (B) for the case of lower particle density in the region. 

A higher value of Tj (run 6) and a lower value of r o (run 7) increase the relative velocity 
of the two shells along with the value of their individual shock propagation Lorentz factors 
(Ff srs ), internal shock energies, magnetic field, 7mi n , and density of higher energy electrons 
injected into the region. The overall flux of the spectrum increases (see Fig. [HD due to the 
reasons cited in (A), (B), and most importantly (C). In both cases, the synchrotron spectrum 
extends into the X-rays: for run 6 v* mQ ~ 6.86 x 10 17 Hz whereas for run 7 the transition from 
synchrotron to SSC component occurs beyond 10 keV. As a result, the spectrum is softer in 
the X-ray range and harder in the Fermi range (see Table HJ) due to reasons given in (D). 
The magnetic field value for these runs is high enough to affect the CDF value according to 
the reason cited in (F). 
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Fig. 4. — Simulated time-integrated SED (left) and lightcurves (right) of our baseline model 
obtained using the input parameters listed in Table [2J The panels on the right show the 
lightcurves that have been evaluated at various energies; from bottom to top, R band, 10 
keV, 1 Mev, and 1 TeV. 




Fig. 5. — Simulated time- integrated SED and lightcurves of a generic blazar source illustrat- 
ing the effects of varying shock parameters, L w (runs 2 & 3), M Q (run 4) according to Table 

El 
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The lightcurves in Fig. [6]have similar profiles as that of run 1, with the same explanation 
as given in the lightcurve profile discussion of run 1 in £0 The pulses of runs 6 (I\ j) and 7 
(r o 1) behave according to the explanation given in (C) at all energies (see Table HJ). 

As shown in Figure [7] and Table [3j increasing the shell widths (run 9) decreases the flux 
and CDF of the spectrum following the explanation given in (G). The values of the peak 
frequencies of both radiation components decrease relative to run 1, thus increasing the SH 
in the X-ray range as explained in (D) but for the case of a low value of 7mi n - On the other 
hand, the location of the trough frequency does not change in this case. As a result, the SH 
for the VHE regime increases as well, as shown in Table [3j 

As can be seen from Figure 0, the lightcurves from run 9 ( Aj & A G f) are much more 
extended than those of run 1. This is because run 9 corresponds to larger widths of the 
colliding shells, thereby resulting in a wider emission region. As a result, the time taken for 
the shocks to propagate through the region is longer. Hence, pulses at all energies follow a 
profile as described in (G). In the case of run 9, the forward shock is the last to break out 
of the system, at t* r ~ 5.27 x 10 5 s, while the reverse shock exits the reverse emission region 
at ~ 5.08 x 10 5 s. 



7.2. Emission Region Parameters 

As can be seen from Fig. [BJ increasing the value of e' e (run 11) increases 7^ in and the 
density of injected electrons in the region. The resultant spectrum and its SH in the X-ray 
and Fermi ranges, along with the value of the CDF (see Table E]), behave according to the 
patterns described in (B) and (D). In the case of decreasing e' e (run 12), the synchrotron 
component peaks in the sub-mm range (^ s * yn ~ 6.21 x 10 11 Hz) with the transition from the 
synchrotron to the SSC component taking place in the optical (f t * urn ~ 6.50 x 10 14 Hz). The 
MeV photons are a result of higher order SSC scattering. The spectrum becomes softer in 
both X-ray and Fermi range (see Table [3]) because the entire spectrum shifts toward low 
frequencies, due to an extremely low value of ■y' min - 

Increasing e' B (run 13) increases the magnetic field while marginally increasing 7^ in 
and the normalization of the injected electron spectrum. The spectral flux, its SH, and the 
corresponding CDF value undergo the same changes as cited in (F), and (D). 

As can be seen from Fig. [HJand Tabled the lightcurve profiles of run 11 (e' e t) follow 
a similar pattern to those of run 1 at all energies. In the case of run 12 (e' e I), since R- 
band photons lie in the highest energy part of the synchrotron component (see Fig. [HJ, the 
high energy electrons responsible for these photons last for a short time in the system. As a 
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Fig. 6. — Simulated time- integrated SED and lightcurves of a generic blazar source illustrat- 
ing the effects of varying shock parameters T; (runs 5 & 6) and r o (runs 7 & 8) according 
to Table [3 




1.0e+06 



v [Hz] obs. time (s) 



Fig. 7. — Simulated time-integrated SED and lightcurves of a generic blazar source illus- 
trating the effects of varying shock parameters Aj & A Q (runs 9 & 10) according to Table 
® 
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result, the corresponding pulse peaks sooner with a smaller FWHM relative to run 1. On the 
other hand, the X-ray photons form the higher energy part of the 1st order SSC component 
and are the result of scattering of high energy synchrotron photons off low energy electrons. 
Thus, the corresponding pulse peaks sooner relative to run 1 but has maximum contribution 
from cooling, and therefore a larger FWHM. The MeV photons form the lower energy part 
of the higher order SSC component and result from the higher order SSC scattering of SSC 
photons off lower energy electrons. Thus, the pulse at MeV energies peaks later and lasts 
for a longer time than that of run 1. The TeV photons are a result of higher order SSC 
scattering off high energy electrons. Since such electrons last for a very short time in the 
system, the pulse has a fairly equal contribution from both acceleration and cooling, making 
it approximately symmetric about its maximum (see Fig. [8]). As can be seen from Fig. [8] 
and Table HI the pulses in run 13 (e' B j) follow the same pattern as explained in (F). 

As shown in Figure [9] and Table [3j increasing C' e (run 15) increases the fraction of 
accelerated electrons in the region while decreasing the value of 7^ in and the density of 
injected high-energy electrons. The effects on the spectrum can be explained using reasons 
given in (E), and (B) for the case of lower density of electrons. As can be seen from the 
figure, the synchrotron component for this case peaks in the IR instead of the optical, as 
in run 1. Thus, the optical photons form the lower energy part of the SSC component and 
are produced from the SSC scattering of much lower energy synchrotron (radio) photons off 
the lowest energy electrons. Hence, the higher order SSC component peaks at ~ 1 MeV 
and is a result of upscattering of first-order SSC photons off low energy electrons. On the 
other hand, the TeV photons are the result of higher order SSC scattering of SSC photons 
off the highest-energy electrons. The spectrum in the X-ray range becomes softer (see Table 
[3]) because, as can be seen from the figure, the transition between the SSC and the higher 
order SSC components lies in the 2-10 keV energy range. On the other hand, the spectrum 
in the Fermi range becomes softer due to the shift in z/g SC to lower frequency. 

As can be seen from Fig. [9] and Tabled increasing the value of a' (runs 17) increases the 
value of 7^, ax while mildly decreasing the value of 7^ in and the normalization of the injected 
electron spectrum. As a result, the SH of the spectrum and the CDF value are hardly 
affected (see Table [3]), but the overall flux of the spectrum increases due to the reasons given 
in (H). 

As can be seen from the right side of Figure |9j pulse profiles of all runs at all energies 
follow the same trend as described in £j7]for run 1. In the case of run 15 ((' e "|), since the 
particle density of the region decreases the t* eak and FWHM of a pulse at a given energy 
follow the same explanation as given in (B) for that case. As can be seen from Fig. [9] and 
Table HI the pulse profiles of run 17 (a' t) are highly similar at all energies to that of run 1, 
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Fig. 8. — Simulated time- integrated SED and lightcurves of a generic blazar source illustrat- 
ing the effects of varying emission region parameters e' e (runs 11 & 12), and e' B (runs 13 & 
14) according to Table [2j 




Fig. 9. — Simulated time- integrated SED and lightcurves of a generic blazar source illustrat- 
ing the effects of varying emission region parameters, (' e (runs 15 & 16), and a' (runs 17 & 
18) according to Table |2j 
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as are their t* eak and FWHM values. 

Increasing the value of q' (run 19) increases 7^ in slightly along with the normalization of 
the injected electron spectrum, while keeping the magnetic field value the same. Thus, the 
overall flux level of the SED rises along with the value of the CDF, as explained in (B). Since 
a softer population of electrons is injected into the region, the peak and trough frequencies 
shift toward lower frequencies, even though the value of 7^ in increases slightly. As a result, 
the SH in the X-ray and Fermi ranges are affected according to (D), but for the case of a 
lower value of 7^ in . 

The lightcurves from run 20 (q' \) (Figure [TO]) show that a harder energy distribution 
of the electron population results in a sharper rise and a steeper decline of the pulse. Since 
the parent population of electrons is more energetic than that of runs 1 (base set) & 19 
(g' f), the electrons lose their energy faster and more energetic synchrotron photons (hard 
X-rays; Figure [TO]) are produced, although only for a short time. As a result, the rise in 
the corresponding pulses at X-ray and optical energies is sooner and the FWHM shorter 
compared to those of runs 1 & 19. For the blazar source under consideration, since TeV 
photons are a result of SSC scattering of higher energy (soft X-rays) photons off higher energy 
electrons, they are produced simultaneously with the lower energy (optical) synchrotron 
photons but with the shortest FWHM (see Table H]) . The TeV photons would lag behind 
the soft X-ray synchrotron photons by a time of the order of the cooling time scale for these 
synchrotron photons. The MeV pulse peaks the last and lasts the longest relative to that of 
run 1 & 19 because for this case, they are a result of SSC scattering of lower energy (optical) 
synchrotron photons off low energy electrons. 



7.3. Jet Parameters 

As can be seen from Fig. [TTJ, increasing R' z (run 21) implies a larger emission region. 
This decreases the overall flux of the SED and the CDF value (see Table [3]) following the 
reasons given in (G). Since the value of 7^ in decreases slightly in this case, the peak and 
trough frequencies, along with the SH, in both energy ranges are affected in the opposite 
way to what is explained in (D). 

For run 21 (R' z f), the pulse profiles at all energies (except at TeV) are much more 
extended relative to run 1 (see Table Hj) due to the reasons explained in (G). Since the width 
of the emission region for both runs 21 & 22 (R' z l) remains the same, the shock crossing time 
does not change. As a result, the TeV pulse lasts only for as long as the shocks are present 
in the system. Thus, the FWHM of the pulse is almost the same as the shock crossing time 
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Fig. 10. — Simulated time- integrated SED and lightcurves of a generic blazar source illus- 
trating the effects of varying emission region parameter q' (runs 19 & 20) according to Table 
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Fig. 11. — Simulated time- integrated SED and lightcurves of a generic blazar source illus- 
trating the effects of varying jet parameter R' z (runs 21 & 22) according to Table [2J 
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(t* r ~ 0.61 days, in the observer's frame) of the forward shock. 

As can be seen from Figure d2] and Table El decreasing the value of 9* hs (run 24) increases 
the overall flux of the spectrum without changing the values of magnetic field, the shocks' 
internal energies, 7mi n > 7ma X > an d the particle density. Since the jet is now closely aligned 
with our line of sight, the overall radiation is boosted more strongly in our direction, thereby 
increasing the overall flux. Due to a lower value of 9* hs , the frequencies when transformed 
into the observer's frame assume higher values. As a result, even though the synchrotron 
and SSC components peak at the same location in the comoving frame, in the observer's 
frame the values of u* yn , z/ t * urn , and z/g SC shift to higher frequencies. This does not affect the 
CDF value, as can be seen from Table |3] for all three runs 23, 24, and 25. However, in the 
case of run 24, this shift of the spectrum towards higher frequencies causes the synchrotron 
component in the observer's frame to extend more into the X-rays, with the transition from 
synchrotron to SSC component taking place at v* urn ~ 4.53 x 10 Hz. As a result, the 
spectrum in the 2-10 keV range appears to be softer as the synchrotron component begins to 
contribute to the soft X-ray flux. On the other hand, the SSC component in the observer's 
frame peaks at v$ sc ~ 1.34 x 10 22 Hz, with the entire hump extending further into the Fermi 
range. This makes the spectrum in the 10 GeV range appear harder, as can be seen from 
the value of ai 0Ge v m Table [3j 

We consider two cases of a larger value of 9* hs . The first one is run 23, where the value 
of #* bs increases such that cos0* bS2g > /3 sh and the line of sight is still within the cone of 
maximum superluminal motion. The second one is run 25 where the value of #* bs is high 
enough to put the line of sight outside the beaming cone. In the comoving frame, this would 
translate into the observer looking at the emission region from behind, as a result of which 
the backside of the region would become visible to the observer before the front side does. 
As shown in Fig. [T2] and Table El for both runs 23 & 25 the opposite trend, compared to 
run 24, is observed. Since the value of # b s in run 25 is higher than that of run 23, the effect 
on the overall flux, the frequency shift, and the SH of the spectrum is even more dramatic 
compared to that of run 23. 

The lightcurves of runs 23 (9* hs t), 24 (9* bs I), and 25 (9* hs f, outside beaming cone) 
are all shifted in time with respect to those of run 1 (see Fig. [T2"j) . In the case of run 24, 
since the value of 9* hs is smaller, the required time for photons at all energies to reach from 
the far side to the near side, of the jet emission region, is reduced. As a result, the pulses at 
all energies peak sooner and last for a shorter duration of time compared to those for run 1 
(see Table H]). Due to a smaller 9* hs , the shock crossing time of the forward shock for run 24, 
in the observer's frame, turns out to be shorter than that for run 1. As a result, the kink 
that is present in the lightcurve of 10 keV photons for run 1, indicating the time of exit of 
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the shock from the system, is smeared out to give a smoother pulse profile at that energy. 
The exact opposite effect can be seen in the pulse profiles at all energies of run 23. In the 
case of run 25, since the line of sight, in the comoving frame of the emission region, passes 
through the region from behind, the light travel time increases by a significant amount. This 
stretches the entire pulse profiles at all energies compared to runs 1, 23, & 24. Once again, 
the value of t* eak for R-band and TeV photons almost marks the time when the shock breaks 
out of the system for all the three runs. 



8. Summary And Conclusions 

We have developed a ID multi-slice internal shock model with radiation feedback scheme 
to simulate the radiative signatures of blazars in a time-dependent manner. We consider a 
cylindrical geometry for the emission region and semi-analytically compute the volume and 
angle- averaged photon escape time scale for it. The inclusion of radiation feedback through 
a multi-slice scheme makes our model a powerful tool to address the issue of inhomogeneity 
in the photon and particle density throughout the emission region. Under this scheme, the 
cylindrical emission region is divided into multiple slices and the subsequent time-dependent 
radiation transfer, within each slice and in between slices, is calculated by using the appro- 
priate photon escape rates of the slices. Our scheme inherently addresses the non-locality 
and time-delay effects due to internal light-travel time while calculating the SSC emission of 
a blazar source. By increasing the spatial resolution along the length of the emitting region 
and considering the contribution of all slices in the emission region, our model mimicks an 
inhomogenous source where the observer sees photons from a mixed population of electrons 
of various ages. Such a scenario is important for cases where the electron population evolves 
on time scales shorter than the dynamical time. Thus our approach facilitates the simulation 
of high-energy variability observed in blazars, especially the SSC dominated TeV sources, 
more accurately. Our model includes synchrotron and SSC emission for evaluating the self- 
consistent time-dependent radiation transfer along the length of the emission region. We 
account for the light travel time delays, in our model, to evaluate the radiative output in the 
observer's frame. 

We consider a single inelastic collision between two shells of plasma of different mass, 
velocity, and energy. The collision results in the formation of a forward and a reverse 
shock, internal to the jet. The relativistic shocks are separated by a contact discontinuity 
(CD). We use the relativistic hydrodynamic Rankine-Hugoniot jump conditions to obtain 
the characteristics of the shock dynamics. The bulk Lorentz factor (BLF) of the emission 
region in the lab frame is evaluated using the equality of pressures of the shocked fluids 
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across the CD. 

The emission region parameters, such as 7 max , 7mim the normalization of the injected 
electron spectrum, and the magnetic field, B, are obtained using the BLFs of the shocks 
and the width of the emission region in the comoving frame of the shocked fluid. The 
emission region parameters are then used to evaluate the radiation spectrum resulting from 
the collision for a cylindrical emission region. A numerical simplification of the synchrotron 
emission formula has been used in order to speed up the computing time. Synchrotron 
spectra calculated with this numerical simplification match those calculated with the exact 
expression to a relative accuracy of better than 1 %. 

The code has been tested to satisfy various analytical constraining conditions. We have 
carried out 25 simulations to perform an extensive parameter study of the code, in order to 
understand the effects of varying various shock (wind luminosity, mass of the outer shell, BLF 
and width of the two shells), emission- region (fraction of shock energy stored in the electron 
and magnetic field energy density, fraction of effectively accelerated electrons, acceleration 
timescale parameter, and particle energy index), and jet parameters (radius of the slice/jet, 
and observing angle inside and outside of the beaming cone) on the resulting SED and 
spectral lightcurves of a generic blazar source. The SED has been parametrized on the basis 
of the change in the level of its overall flux, shift in its peak and trough frequencies, change 
in its spectral hardness in the 2-10 keV and 10 GeV energy ranges, and the Compton 
dominance of the SSC component over the synchrotron component. The lightcurves at 
optical, X-ray, MeV, and TeV energy ranges have been parametrized on the basis of the time 
of exit of the shock from the emission region, the time of the peak, and the FWHM of a 
pulse at a given energy. 

As demonstrated in §HJ both SEDs and lightcurves give us an idea about the dominant 
radiation mechanism responsible for producing radiation of a particular frequency and the 
corresponding time lags between various frequency bands. As shown in §|6l our model can be 
exploited to reproduce long and short-term variability of blazars based on the combination of 
shock, emission-region, and jet parameters. The symmetric and asymmetric profiles observed 
in the lightcurves of blazars are frequency dependent. Lightcurves at higher synchrotron and 
SSC frequencies tend to be more symmetric than th eir lower frequency counterparts (see 



e.g., iChiaberge fc Ghisellinilll999l ; ISokolov et all 120041 ) for the reasons explained in £j6j As 
indicated in §|6j such profiles of the lightcurves can be reproduced using the right combination 
of shock, emission-region, and jet parameters. The symmetric and asymmetric profiles of the 
lightcurves are also dependent on the angle at which the jet is being viewed (see Fig. [T2"|) . As 
mentioned below, since we are not carrying out a 2-D simulation of the jet radiation processes, 
we are limited by the internal light crossing time in the radial direction. As a result, the 
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symmetric profiles that appear from viewing the jet at non-zero angles ( ISokolov et al.ll2004[ ) 
and are dominated by light crossing time effect across the emission region are not accounted 
for. 

The particle injection spectra can be used to identify the characteristic signatures of var- 
ious particle acceleration mechanisms. The particle acceleration at parallel shocks can give 
rise to electron spectra following a power-law dis tribution with an injection index of q' ~ 2.2 - 
2.3 (lAchterberg et al.ll200ll ; iGallant et al.lll999l ). Oblique shocks, on t he other hand, are be- 
lieved to produce much softer inj ection spectral indices with q' > 2.3 (jOstrowski fc Bednarz 
20021 ; iNiemiec &: Ostrowskil 120041 ) . In contrast, 2nd order Fermi acceleration might result in 
a harder injection index of the order of q' ~ 1 or beyond if the plasma is turbulent, behind 
the shock front. Thus, the value of q' used for modeling the SED and spectral lightcurves of 
a blazar source can give us insights on the kind of shock propagating through its jet and the 
dominant mode of acceleration present in the source. Our future modeling efforts will help 
us gain a deeper understanding of the real data. 

We would like to point out that for simplicity we have considered the role of the magnetic 
field to be negligible in the dynamics of the shock energetics. The structure of the mean large 
scale magnetic field plays an important part in particle acceleration and shock dynamics. 
We also emphasize that for making the problem at hand more tractable we chose to have 
the particle injection spectral index, free parameter. The value of q should really 

be the result of shock dynamics and the underlying magnetic field structure. In our model, 
the escape timescale for electrons from the emitting region is a weakly constrained, free 
parameter. It involves a fudge factor that parametrizes the randomness of the magnetic field 
direction. Hence a more detailed model of the radiation transfer in blazar jets, dealing with 
particle accleration, must address these physical parameters. 

We note that our model, in its current form, does not address the internal light travel 
time in the radial direction. The time resolution of our model is limited by the light crossing 
time across the cylinder. This feature is more relevant for carrying out a comparison between 
synchrotron and SSC pe ak times at highest frequencies, especially for the case where the jet 
is being viewed face-on (jSokolov et al.ll2004[ ). A 2-D multi-celled approach would be a more 
appropriate tool to incorporate this aspect of frequency-dependent time lags. 
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A. 



Coefficients of T. 



The coefficients of Eqn. [TT] are as follows: 



a = 7 2 [(Y - l) 2 + 4Y (n - 2x 2 m)] , 
b = 2^g [(Y + 1) (r o + YTi) - 2Yxpm] , 



c = 27Y (2 - 37) n + 27 (7 - 2) r - 2Yxm (# 2 - 47 2 x) - 4Y^g + (l - 7 2 ) (l + Y 2 ) , 
d = 2(7 [Yp (1 + 7 (x (1 - 2q) - 1)) - (1 + 7) (r + Y 2 r rmj ) + 7 (rH + Y 2 rf)] , 
e = 7 2 (r^ + Y 2 rf) + 27Y (2 - #72 - 7 (1 + x 2 )) + r (l - 7 2 ) - 2g 2 Yxq - 2Y . (Al) 

Here, x = V Q Y h m = 1 - f3 Q f3 h q = f3 Q f3 h n = T 2 + T 2 , p = T + T h r = T 2 Q + Y 2 T 2 , 



g — 1 — 7, and Y = p-J p Q . All unprimed quantities refer to the AGN frame and quantities 
with an overline refer to the unshocked fluid frame. 



All quantities mentioned here refer to the comoving frame of the emission region. The 
intermediate steps used to derive equation fl27j) are the following: 

As shown in Figure [TBI (a), / CT it+ is the distance traveled by a photon that escaped from 
the edge of the region in the forward direction, at an angle # C rit+ with the z-axis, and similarly 
for / crit _ at an angle n — 9 crit - with the z-axis; b is the corresponding projection of these 
lengths on the horizontal plane (Figure [TBI (b)). 

Using the law of cosines and solving the quadratic equation in b, and realizing that 
negative distances are unphysical, an expression for b can be obtained as 



B. Photon Escape Time Calculation 




(Bl) 



Now, referring to Figure [TBI we can write 



tan 6, 



b 



, tan 9, 



-b 



crit+ 



h — Z 



crit— ) 

Z 
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h sc - = — for -1 < \i < /i crit -, 
//_ 

; b r 

< psr — , IOr jJ>crit- < fi < /X C rit+> 



"esc-f 



h — z 



for // crit+ <//<+!, 



/^crit+ 



h — Z 



b 2 + (h- zf 



A^crit- 



(B2) 



Substituting equations (1B2j) in equation ( )26|) and then equation ( )26|) in (|25|) . equation 
( |25|) can be written as 



2tt 



(tph ' esc ' v) " ( ^ )( 4^ } 

Merit— Mcrit+ +1 

1 Merit- /^crit+ 



2iv R h 








• rdr dz 



which gives 



(B3) 



(^ph,esc,v) 
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>irR 2 h JK 4nc J 



2n R h 
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Mcrit- 



A t crit + 



+ 1 



— -d[i- 



/^crit - 



'1 



-.dfj, + 



h — z 



djj,- 



/XcritH 







i rdr dz 
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Carrying out integration over /i, and realizing that the integration of the right hand side 
expression over the angle $ is going to give 2 ir, the above equation can be written as 



R h 



(Wsc,v) = J^h J {J(^ z ~ h ^ - z) - z\nz}+ h + I 2 + I 3 )dzjrdr , (B5) 







where 
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2- 



/, = I \ n (b 2 + (h - z) 2 )d(f>, 



2tt 



57TC 



ln(6 2 + z z )d<p, 



271 



^6[arcsin( — ) + amsin( )\d<j> 



4nc 



^b 2 + (h- z)< 



Vb 2 + z 2 ' 



(B6) 



The above three integrals can be solved semi- analytically to obtain the final expression 
for the volume and angle-averaged photon escape timescale for a cylindrical region given in 
Eqn. d23 



C. Synchrotron Function Fitting 

All quantities mentioned here refer to the comoving frame of the emission region. The 
function R(x), used in the ca lculation of the synchrotron photon production density rate 
(see §|H Eqn. [51]) . is given by ( ICrusius fc Schlickeiserlll986f ) 



R(x) 



7TX 



(CI) 



where x 



A-Km e c v 
3cB" , < 



is the nor malized frequency, and W\ tfl (x) denotes Whittaker's function 



( lAbramowitz fc Stegunlll970l ). 



In order to save CPU time while computing equation ( 13 ip O, we numerically approx- 
imate the function, R(x), as follows: 



R(x) = C lX Pl e- x - C 2 x~ P2 e- x , 
where C x = 1.08895, C 2 = 2.35861 x 10~ 3 , pi = 0.20949, and p 2 = 0.79051. 



(C2) 



The asymptotic expansions of Whittaker functions fjCrusius fe Schlickeiserl Il986f ) that 
are used in the calculation of Eqn. |3U for small (i C 1) and large (x ^> 1) values of x, are 
no longer required when using the above approximation. Since R(x) has been approximated 
for the expression given in Eqn. ( ICip . we normalize R(x) at lower values of x, where R(x) 
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oc x 1 / 3 . This is carried out to obtain a good fit to the exact expression for the entire range 
of x values. 

Figure O shows the comparison of simulated SEDs of a generic blazar source, obtained 
from using the approximation and the exact expression. As shown in the figure, the approx- 
imation is accurate to better than 1% over the entire frequency range. 
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Fig. 12. — Simulated time-integrated SED and lightcurves of a generic blazar source il- 
lustrating the effects of varying jet parameter 9* hs (runs 23, 24 & 25) according to Table 

m 




Fig. 13. — In (a) we illustrate the three angles of escape along with their corresponding 
escape lengths. In (b) we show the angle consideration and the projection of the escape 
lengths on the horizontal plane of the cylinder, when viewed from the top. 
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Fig. 14. — Comparison of the exact expression and numerical approximation of R(x). The 
fit is accurate to ~ 0.5% for the entire frequency range. 



